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ABSTRACT 



Aims. The aim of this paper is to study the vertical profile of small dust particles in protoplanetary discs in which angular momentum 
transport is due to MHD turbulence driven by the magnetorotational instability We consider particle sizes that range from approxi- 
mately 1 micron up to a few millimeters. 

Methods. We use a grid-based MHD code to perform global two-fluid simulations of turbulent protoplanetary discs which contain 
dust grains of various sizes. 

Results. In quasi-steady state, the gravitational settling of dust particles is balanced by turbulent diff'usion. Simple and standard mod- 
els of this process fail to describe accurately the vertical profile of the dust density. The disagreement is larger for small dust particles 
(of a few microns in size), especially in the disc upper layers (Z > 3H, where H is the scale-height). Here there can be orders of 
magnitude in the disagreement between the simple model predictions and the simulation results. This is because MHD turbulence is 
not homogeneous in accretion discs, since velocity fluctuations increase significantly in the disc upper layer where a strongly mag- 
netized corona develops. We provide an alternative model that gives a better fit to the simulations. In this model, dust particles are 
diff'used away from the midplane by MHD turbulence, but the diff'usion coefficient varies vertically and is everywhere proportional to 
the square of the local turbulent vertical velocity fluctuations. 

Conclusions. The spatial distribution of dust particles can be used to trace the properties of MHD turbulence in protoplanetary discs, 
such as the amplitude of the velocity fluctuations. In the future, detailed and direct comparison between numerical simulations and 
observations should prove a useful tool for constraining the properties of turbulence in protoplanetary discs. 

Key words. Accretion, accretion discs - MHD - Methods: numerical 



1. Introduction 

Thanks to the Spitzer Space Telescope, the last few years have 
seen a dramatic improvement in our knowledge of dust emission 
I features arising at mid-infrared wavelengths from protoplanetary 
' discs surrounding brown dwarfs, T Tauri stars and Herbig Ae/Be 
stars. The properties (size, composition) of these grains can be 
studied in great detail by analyzing the shape and strength of 
' these emission features. The main result is that dust particles 
have been significantly processed compared to their interstellar 
medium cousins: grains are bigger (up to a few microns) and, 
for reasons not yet fully understood, crystallinity ap pears to be 
common among all observed spectral types ( A pai et alJ 120051: 
iKessler-Silacci et alJ 120061: iFurlan et al]l2006h . Because proto- 
planetary discs are optically thick, this mid-i nfrared emission 
mostl y arises from the upper layers of the disc (Ivan Boekel et al.l 
'2003VDullemond & Dominik"2004' '2008), in the so-called "su- 
perheated" layer of Chiang & Goldreich ( 1997), and is attributed 
to the disc inner radii. In the case of T Tauri stars of solar-type, 
for example, the dust emission zone lies within a few tenth s of an 
astronomical unit (AU) of the central star (Kessler-Silacci et aP 
I2Q07I) . 

The presence of solid particles high above the disc mid- 
plane is the signature of the turbulent nature of the flow in 
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protoplanetary discs. Turbulent velocity fluctuations lift up dust 
particles that would otherwise settle toward the di sc midplane 
(iDubrulle et aDll995l: IPullemond & Dominikll2004l) because of 
the vertical component of the central star's gravitational po- 
tential. The vertical profile of the dust density is thus deter- 
mined by the balance between gravitational settling toward the 
equatorial plane and upward lift due to turbulence. Its typi- 
cal scale-height is a function of grain size. Small particles are 
well coupled to the gas, and can be transported to higher al- 
titudes (or smaller gas densities) above the equatorial plane 
than larger particles. Although th e issue is still under discus- 
sion (iDuUemond & Dominikl2008b . this difl'erential settling pro- 
cess can be inferred from a detailed analysis of the observa- 
tions. This is the case for observations carried out by Spitzer 
alone (Furlanetal. 2005; Sicilia-Aguilar et al."2007'), by com- 
bining Spitzer observati ons with observations at other wave- 
lengths (iPinte et al.ll2008h or even using completely differe nt ob- 
servational strategiesT Pinte et al.ll2007l:lRettig et al.ll200 6l). Dust 
gravitational settling is also known to alfect the spectral energy 
distribution of protoplanetary di scs, as shown for exa mple by 
[DuUemond & Dominik (2004) or lD' Alessio et al.1 (12006) . 

In general, such observational diagnostics of gravitational 
settling are inferred using simple parametric presc ription s of the 
effect of turbulence or of its consequence. Pinte et SI (|2007i 
I2008h assume that the vertical profile of the dust density is a 
Gaussian (as is the case for the gas in the isothermal limit), 
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leaving open the possibility that the dust disc scale-height de- 
pends on the size of the particles. The point of their analysis is to 
demonstrate that this is indeed the case in GG Tau and IM Lupi. 
iFurlan et al.l (l2QQ5l) use the models developed by D' Alessio et ajj 
and assume a constant dust-to-gas ratio for two popula- 
tions of solids, one consisting of small and well mixed parti- 
cles, the other composed of large particles. Using Spitzer spec- 
tra, they demonstrate that the latter is depleted in the disc upper 
layers for the vast majority of a sample of 25 stars in Taurus. 
iDullemond & DominikI (120041) use a more physical approach 
in which turbulence is modeled as a diffusive process with a 
spatially constant diffusion coefficient (see also Dubrulle et al. 
P^5; Schrapler & Henning 2004). Using this formulation, one 
can derive analytic expressions for the dust density vertic al pro- 
file. This is also the approach followed by Rettig et^aP (l2006l) 
who used the strong settling (thin dust disc) limit of these for- 
mulae. 

All of these approaches are very useful since they provide 
analytical formulae or a small number of model parameters to 
fit. They can thus be incorporated into radiative transfer tools, 
which can then generate a large grid of models to be used for 
interpreting the observations. All of them establish, at least qual- 
itatively if not quantitatively, the basic result that differential 
gravitational settling occurs in protoplanetary discs. However, 
none of these derivations follows directly from first principles, 
but instead rely at best on a set of unchecked assumptions con- 
cerning the properties of the underlying turbulence. The origin 
of the latter, however, is thought to be magnetohydrodynamic 
in nature and driven by t he ma gnetorotational instability (MRI, 
iBalbus & Hawlevl ll99l[ [1998). It is now possible to perform 
global numerical simulations of turbulent protoplanetary discs, 
including solid particles, and to use these simulations as a test of 
the models described above. This is the purpose of the present 
paper. Despite the recent increase in computational resources, 
it is important to keep in mind that such simulations are still 
extremely challenging and need to use as simplified a set-up 
as possible. This is not without consequences for their realism. 
Important issues related to the MRI, such as small scale dissi- 
pation (Fromang et al. 2007; Lesur & Lonsaretti 2007) and the 
possible presence of a dead zone (Gam mielll996>) still have to be 
ignored. We will return to these issues in the last section of the 
paper when discussing the limits of our work. 

Of course, the effect of MRI-induced MHD turbulence on 
dust dynamics and dust vertical settling itself has already been 
studied in numerical simulations, but this is the first time that 
both effects are incorporated in a singl e global simulation that 
takes vertical stratification into account. iBarriere-Fouchet et al.l 
(2005) performed global simulations of dust settling but ne- 
glected the effect of MHD turbulence. This prevents the sys- 
tem from reaching a quasi-steady state. Other global simu- 
lations, taking MRI-driven turbulence into account and in- 
cluding dust particles, ne glected gaseous vertical stratification 
(iFromang & Nelson'2005') and instead focussed on the radial mi- 
gration of larger bodies. Ly ra et al.„ (2008) also neglected the ver- 
tical component of gravity acting on the gas, but included its ef- 
fects on the dust particles. While this approach produces settling 
of the dust particles toward the midplane, its neglects the spatial 
inhomogeneity of the turbulence induced by vertical stratifica- 
tion of the gas. We shall see in the course of this paper that the 
latter is important when considering dust settling of small parti- 
cles. The other published numerical simulations studying the ef- 
fect of turbulence on dust dynamics were local simulations that 
use the shearing box model ( Goldreich & Lynden-Belll 11965). 
A large number of them neglected density stratification in the 



vertical direction (jJohansen et al. I l2006al: Ic^alfido et al.ll2005l: 
Johansen et al. 1 12006b') and thus only studied dust diffusion in 
an hom ogeneous en vironment. Others (Fromang & PapaloizoiJ 
I2OO6I: ICarballido et al.l t2006) included vertical density stratifica- 
tion but considered particles larger than one millimeter. Here, 
we want to concentrate on smaller particles that produce an ob- 
servational signature at mid-infrared wavelengths. We note that 
such a simulation could in principle be done in the framework 
of the shearing bo x model. Indeed, this was done recently by 
iBalsara et al.l (12008 ). although the vertical extent of their shear- 
ing box is smaller than the simulations we present in this paper 
and thus less appropriate to study the dust distribution in the disc 
corona (i.e. above three scale-heights). However, since this work 
is intended to mark the beginning of an effort to compare obser- 
vations and numerical simulations directly, it makes more sense 
to compute global models as these will be more readily compa- 
rable with the observations as they become more realistic. Our 
strategy in this pap er will be simple. We wil l use exactly the 
set-up presented by lFromang & NelsonI (l2QQ6l) for global simu- 
lations of turbulent and stratified protoplanetary discs. Dust par- 
ticles of various sizes will be added to the disc and their sub- 
sequent evolution will be analyzed and compared with simple 
models of dust stratification in protoplanetary discs. 

The plan of the paper is as follows. In Sect.O we introduce 
the model we use for the disc as well as convenient dimension- 
less parameters that appear in the problem. The relationship be- 
tween these quantities and physical parameters in real systems 
will be outlined. In Sect. O we describe more quantitatively the 
different ways to model the quasi-steady state dust distribution 
resulting from the balance between dust settling and turbulent 
diffusion. These models are then compared with the results of 
our numerical simulations in Sect.lH Finally, in Sect. [5] we dis- 
cuss the implications and limitations of our work, and point the 
way toward future improvements. 

2. Definitions 

In this section, we describe the general properties of the disc 
model and the dust parameters we used. We also introduce the 
mathematical notation that will be required in the following sec- 
tions. 

2.1. Coordinate systems 

In this paper, we will use both cylindrical and spherical coordi- 
nate systems. The former will be denoted by (7?, 0, Z), and will 
be used mostly in the present Sect. [2] and in Sect. [3l Spherical 
coordinates will be used when we describe the numerical simu- 
lations and will use the notation (r, 6, cp). 

2.2. Disc model 

We consider a disc extending in radius between an inner radius 
Rdjn and an outer radius Rd,out- For simplicity and computational 
reasons we define the initial disc structure using straight-forward 
analytic functions. The equation of state is locally isothermal: 
the sound speed, c^, only depends on R and is constant in time. 
Both Cs and the disc midplane gas density, p^, obey power laws 

c^siR) = , (1) 

/ ^ \-3/2 

PmidiR) = PO — , (2) 
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where cq and po are the sound speed and the midplane gas den- 
sity at a radius Rq, respectively. The disc is initially axisymmet- 
ric and in radial and vertical hydrostatic equilibrium. The spatial 
distribution of density p(R, Z) and angular velocity Z) can 
thus be approximated by 

p(R, Z) = pMe-''^'"' = PO I e-^'^'^' , (3) 



-3/2 



(4) 



in which Qq = ^GM/7?q. The disc scale-height, H, is given by 

(5) 

where Hq = cq/Qq is the disc scale-height at Rq. 



2.3. Integrated quantities 

The density distribution can be used to work out the surface den- 
sity of the disc: 



(6) 



where the second relation serves as a definition for So- The total 
disc mass follows by radially integrating the surface density 
between R^jn and Rd,out' 



471^ 

T 



^0 



3/2 



(7) 



where we have assumed Rd,m ^ Rd,out- This relation can be used 
to express the disc surface density at Rq as a function of the disc 
parameters and R^: 



^0 = 



3M^ / 7?o 



AnRl \Rd,out 



3/2 



(8) 



2.4. Dust size 



In this paper, we shall study the eff'ect of MHD turbulence on the 
dust. Gas aff'ects solid body dynamics through the drag force it 
exerts on the dust particles. In the Epstein regime we are inter- 
ested in, this force Fa takes the simple form 



Fd = - 



V-Vd 



(9) 



where v and Vd are the gas and dust velocities, respectively, is 
the dust stopping time. This is the typical time it takes for dust 
particles initially at rest to reach the local gas velocity. It depends 
on the dust particle mass density ps and its size a through 



Psa 

PCs 



(10) 



A relevant dimensionless parameter in the problem is the quan- 
tity flTs. It compares the stopping time to the dynamical time. 
When flTs ^ 1 , the stopping time is much smaller than the or- 
bital period Torb and the dust essentially follows the gas. When 
Q.Ts ~ 1 or larger, becomes compara ble to Tprh and dust and 
gas sta rt to decouple. As pointed out by lDullemond & DominikI 
(I2QQ4 . this occurs at all radii for a given particle size a provided 



Z is large enough. Using the disc parameters introduced above, 
can be expressed as a function of position according to 



poCo\Ro) \Ro) 

where the parameter (^Ts)o is the value of Qr^ aiR = Rq in the 
disc midplane. It can be expressed in terms of the disc surface 
density at Rq using Eq. ©i 



1/2 



(Qt,)o= ^^2^^. 

^0 



(12) 



Using this expression along with Eq. ([8]), it is possible to express 
the dust size in term of (^Ts)o, the disc parameters (mass and 
radius) and Rq: 



3Md iRd 



2^fhlpsRl \ ^0 



-3/2 



(^^T,)o. 



(13) 



2.5. Converting to physical units 

Eq. ([5]) and ([T3]) can be used to convert the dimensionless quanti- 
ties describing the problem into physical quantities. When doing 
so, the numerical simulations we will describe in Sect. IH should 
be thought of as covering a small fraction of the total disc, going 
from Rin = Ro> Rd,in to an outer radius Rout < Rd,out- 

The dimensionless parameters describing the disc and dust 
particles, and the results of the numerical simulations presented 
below, can be reseated to any physical system upon specifying 
the disc mass, the outer radius of the disc, and the value of R^ 
(in astronomical units). For example, the disc surface density 
can be written as 



Md 



O.OlMc 



Rd,out \ 

300 Au] 



-3/2 



(^1 

\IAUI 



-1/2 



gem 



(14) 



Likewise, taking p^ = 1 gem ^ and using Eq. ([T3l) , we obtain an 
expression for the dust size: 



163 



0.01 /\0.01M( 



Rd,out \ 

300 AU I 



-3/2 



^0 y 

lAu) 



■1/2 



yum. (15) 



In Sect, m we will describe the results of three simulations. 
They are characterized by (^Ts)o = 10"^, 10"^ and 10""^. In a 
disc of mass 0.0 IM© and 300 AU in size, they would correspond 
to dust particles of size 163, 16 and 1.6 yum, respectively, if we 
take Ro = 1 AU. (This would mean that we consider the simu- 
lation to cover the radial extent 1 to 8 AU.) For the same disc 
mass and size, if we now take Ro = OA AU (i.e. we consider the 
simulation to cover radii ranging from 0.1 to 0.8 AU), the three 
sizes are 500, 50 and 5 jim. Note, however, that these numbers 
are only illustrative. In general, for a given value of Rq, the size 
of the dust particles decreases when the disc mass is decreased 
or its outer radius is increased. 



3. Dust settling in turbulent discs 

As pointed out in the introduction, a steady state is reached in 
a turbulent protoplanetary disc in which turbulent fluctuations 
oppose and balance against dust settling. In this section, we de- 
scribe three approaches that can be used to model the vertical 
profile of the dust density. 
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3.1. A Gaussian profile 

The simplest approach, and one tha t is commonly u sed when 
attempting to interpret observations ("Pinte e t aLll2QQ8 h, is to as- 
sume that the dust density follows a Gaussian distribution, as the 
gas does, but with a different vertical scale-height Hd'. 



Pd - Pd,mide 



(16) 



where pd,mid is the dust midplane density. In this approach, Hd is 
different for each dust particle size a. For example, while trying 
to model the dust pro perties in the prot oplanetary disc orbiting 
the M dwarf IM Lupi Jpinte et~aD (l2QQ8l) found Hd oc a"^ One 
purpose of this paper is to establish whether such a description is 
supported by numerical simulations of turbulent protoplanetary 
discs. 

3.2. Turbulence as a diffusive process 

A more physical approach is to describe the transport of 
the dust particles by the turbulent fluctuations as a diffu- 
sion process. This has been u sed commonly in the literature 
(iDullemond & DominikI l2004l) . In this approach, the vertical 
evolution of the dust density ca n be described by the fol- 
lowing partial differe ntial equation (ISchrapler & Hennili^l2QQ4l: 
iDubrulle et aT]ll995l) : 



dpd ^ / o2 N ^ 
— - — (ZQ Tspd) = — 



dt dz 



dz 



dz\p 



(17) 



where pd is the dust particle density and Z) is a diffusion co- 
efficient that quantifies the turbulent diffusivity. This equation 
models the balance between vertical settling and turbulent dif- 
fusion. When looking for a steady state vertical profile for the 
density, the time derivative vanishes. Upon integrating once and 
rearranging terms, Eq. (fTTl) gives 



D 



-z. 



(18) 



The vertical integration of the last equation requires the knowl- 
edge of the diffusion coefficient as a function of z. The simplest 
solution is to assume that it is constant. This is the approach 
described in the following subsection, while in Sect. 13.2.21 we 
outline a possible alternative. 

3.2.1 . A constant diffusion coefficient 

When the dust diffusion coefficient D is constant, Eq. ([TS]) can 
be integrated to give 

2/f 2 ) ) 2//2 



Pd = Pd,mid exp 



D 



mid 



exp 



(19) 



where D is a dimensionless diffusion coefficient defined as Z) = 
bcsH. The quantities pd,mid and (SlTs)mid only depend on R and 
are to be evaluated in the disc midplane. Note that, while deriv- 
ing the last equation, we have assumed that the vertical distribu- 
tion of gas remains Gaussian at all times, in agreement with local 
shear ing box numerical simulations of the MRI ( Miller & Stone 
[200QI) . 

A common practice in this context is to exp ress D in units 
of the standard a parameter introduced by Sha kura & Sunvaev 
D is the n written as follows (iDullemond & DominikI 
2004; Sc hrapler & Henning..2004) : 



Sc 



(20) 



where Sc is the Schmidt number. In zero net flux MHD 
turbulence, the Schmidt number has been measured to 
be of order unity in local simulations of unstratified 
(iJohansen & Klahr 2005; Johansen et al. 2006b) or stratified 
discs (iFromang & Papaloizou 2006: Jlgner & Nelson .20081 In 
the present paper, we will tune the Schmidt number in order to 
obtain the best agreement with the numerical simulations. 



3.2.2. A vertically varying diffusion coefficient 

Dust particles are diffused away from the disc midplane by the 
turbulent velocity fluctuations. Thus, the dust diffusion coeffi- 
cient is intimately linked to the turbulence properties and partic- 
ularly to the gas velocity fluctuations. It would then seem natural 
for D to be constant in space if the turbulence was homogeneous. 
However, because of the vertical stratification, MHD turbulence 
is not homogeneous in protoplanetary discs and it is very likely 
that D varies with Z at a given radius (even in the absence of 
a dea d zone). Using local vertically stratified simulations of the 
MRI, Fromang & Papaloizou| (l2006l) showed that the following 
simple expression gives a fairly good estimate to the numerically 
derived diffusion coefficient: 



(21) 



In this equation, Svz stands for the turbulent velocity fluctu- 
ations and Tcorr is the correlation time of these fluctuations. 
Numerical estimates of both terms thus provide a path to cal- 
culating the value of D. Their vertical variations will be investi- 
gated in Sect. HI 

Of course, the drawback of this approach is that it becomes 
impossible to explicitly integrated Eq. (fTSl) . We cannot provide 
an analytical expression for the vertical profile of the dust den- 
sity and shall rely on a numerical integration once the vertical 
profile for D is extracted from the numerical simulations. 



4. Numerical simulations 

4.1. Set-up 

The simul ations presented in this paper are run using the code 
GLOBAL (lHawlev&Stonelll995h . which solves the ideal MHD 
equations using a spherical coordinate system as defined in 
Sect. 12. 1| The set-up we used is exactly that of model S5 in 
iFromang & NelsonI (l2006l) . Here, we describe it only briefly. 

At the start of the simulation, the disc model presented 
in Sect. 12.21 is initialized on the grid. The units are such that 
GM = 1, Co = 0.1 andpo = 1 (i.e. H/R = 0.1 throughout the 
disc). The computational domain covers the range 7?/^ — Rq — 1 
to Rout = 8 in radius and the interval [0,7r/4] in 0. In the 6- 
direction, the grid extends to 4.3 scale-heights on both sides 
of the disc equatorial plane. The resolution for each simulation 
is {Nr^N^.Ne) = (364, 124. 213). Following From ang & NelsonI 
( 2006), a weak toroidal magnetic field is added to the disc, and 
small random velocity perturbations are also imposed. Time is 
measured in units of the orbital period at the inner edge of the 
computational domain in the following sections. 

Because of the MRI, the presence of a weak magnetic field, 
together with the velocity perturbations, begins to drive MHD 
turbulence and angular momentum transport through the disc 
within a few orbits of the simulations starting. To reach a mean- 
ingful quasi steady state, however, the model is first evolved for 
430 orbits without dust particles. At that stage, the gas density 
is reset to its initial distribution and dust particles are introduced 
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Fig. 1. Snapshots of the gas density (left panel) and the dust density {right panel) aXt = 600 for the case (^t)o = 0.01. 



such that the dust-to-gas ratio is uniform through the compu- 
tational domain (note that we neglect the back reaction of the 
solids onto the gas, so that the value of this ratio has no phys- 
ical consequences). We ran three simulations for three different 
particle sizes. The three values of (^Ts)o associated with these 
sizes are 10"^, 10"^ and 10"^. All the simulations are further 
integrated for about 200 orbits until the dust distribution itself 
reaches a steady state in which gravitational settling is balanced 
by turbulent diffusion. Examples of the disc structure at the end 
of such a run are illustrated in Fig. [T] Two snapshots of the gas 
(left panel) and dust (right panel) density in the (R, Z) plane are 
shown in the case (Qt^)o = 0.01 at time t = 600. The dust disc 
appears thinner than the gas disc, indicating that significant set- 
tling has occurred. 

In the following subsections, we will compare the dust dis- 
tributions we obtained for the different sizes to the models de- 
scribed in Sect. [3l To do so, all relevant physical quantities will 
be averaged in time between ^ = 550 and t = 600 using 50 snap- 
shots so that they become statistically significant. We first start 
by describing the relevant properties of the turbulence that re- 
sult from this procedure before concentrating on the degree of 
settling as a function of size. 
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Fig. 2. Radial profile of a (solid line), a Max (dashed line) and 
aRey (dotted line) for the case (^t)o = 0.001. The data have 
been averaged in time between ^ = 550 and t = 600. 



4.2. Turbulence properties 

The first relevant property of the turbulence is the vertically av- 
eraged angular momentum transport it gener ates. It is commonly 
measured using the parameter a. Following Fromang & Nelson! 
(i2006) , we calculated a as a function of radius according to 



a = aRey + a Max = ^ = , (22) 

P 

where aRey and a Max correspond respectively to the Reynolds 
and Maxwell stress contributions to a. The overbar symbols 
denote density- weighted azimutha l and vertical averages (see 
Eq. (6) of Fromang & Nelsonll2006 ). To reduce statistical noise, 
a was further averaged in time between ^ = 550 and t = 600. Its 
radial profile is shown in Fig. [2] for the model having (Q.t)o = 
0.001. The dashed and dotted lines, respectively, show the vari- 



ations of a Max and aRey, while t he solid line represents a, the 
sum of the two. As described in iFromang & NelsonI (l2006h . a 
presents large oscillations in space and time, but its time aver- 
aged value is well behaved, varying between ~ 5 x 10"^ and 
~ 1 .5 X 10"^, in agreement with the results of Fromang & NelsonI 
( 2006). As is usually obtained in numerical simulations of this 
type, the Maxwell stress dominates over the Reynolds stress by 
a factor of about two to four. The volume averaged value of a is 
also in agreement with the results of Fromang & Nelson ( 200^. 
Indeed, we obtained = 5.6 x 10"^ 5.5 x 10"^ 7.5 x 10"^ 
7.1 X 10"^ 6.4x 10"^ and 6.6x 10"^ respectively at times t = 550, 
560, 570, 580, 590 and 600. For the duration of the simulation, 
the turbulence is clearly in a quasi steady state. 

In this paper, we present three simulations for different parti- 
cle sizes. Because these simulations were obtained under differ- 
ent computing set-ups (i.e. using different numbers of CPUs), the 
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Fig. 3. Radial profile of a (time averaged between ^ = 550 and 
t = 600) for the models (Qt)o - 0.01 (solid line), 0.001 (dotted 
line) and 0.0001 (dashed line). Angular momentum transport is 
similar at all radii in the three models, despite differences arising 
because of the stochastic nature of MHD turbulence (see text for 
details). 



0.4 




Fig. 4. Vertical profile of the vertical velocity fluctuations (in 
units of the speed of sound) for the model having (^t)o = 0.01 
(solid line), 0.001 (dotted line) and 0.0001 (dashed line) at 
R = 2.93. The simulations data have been averaged between 
t = 550 and t = 600. 



details of the turbulent flows differ from one run to another. This 
is s imply due to the chaotic nature of turbulence ( Winters et al. 
12003 ). However, the statistical properties of the turbulence are 
similar. This is shown on Fig. [3] where we compare the radial 
profile of a for the models (^t)o = 0.01 (solid line), 0.001 (dot- 
ted line) and 0.0001 (dashed line). For each case, the curves are 
averaged in time between ^ = 550 and t = 600. Fig. [3] demon- 
strates that we obtain very similar values of a in the different 
simulations. Thus it is meaningful to compare the dust distribu- 
tions in the three different models. 

As explained in Sect. 13.2.21 the vertical velocity fluctuations 
and the correlation timescale of the turbulence are also of im- 
portance in affecting the diffusion of solid particles. In Fig. IH 
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Fig. 5. Time variation of the correlation function of the vertical 
velocity fluctuations in the disc midplane (solid line) and in the 
disc corona (dashed line). The dotted lines represent functions 
decreasing exponentially toward zero with typical times tq = 
0.05,0.1, 0.15 and 0.2 orbit. They can be used to estimate the 
typical correlation timescale of the turbulence. At both locations, 
it is not far from the value Tcorr = 0.15 orbit we used when 
modeling the results of the numerical simulations. 



we show the vertical profile of the vertical velocity fluctuations 
at 7? = 2.93, normalized by the speed of sound and averaged in 
time between ^ = 550 and t = 600 for the models (I1t)o = 0.01 
(solid line), 0.001 (dotted line) and 0.0001 (dashed line). Again 
the results are very similar for each model and in agreement with 
those of Fromang & NelsonI (l2006h : the average Sv^ is of order 
5% of the speed of sound in the disc midplane before rising up 
to values of the order of 20 - 30% in the disc corona, where 
weak shocks develop in l ocations where convergent flow speeds 
exceed the sound speed (iFromang & NelsonI l2006l) . These su- 
personic turbulent motions are driven by magnetic stresses in 
regions where the Alfven speed exceeds the sound speed. Fig. [4] 
shows that Sv^ varies by a factor of about 5 between the disc 
midplane and its corona. 

As expressed by Eq. (|2T1) , the dust diffusion coefficient also 
depends on the turbulence correlation timescale Tcorr- Although 
D depends on the value of Tcorr to the first power only, while it 
depends on the velocity fluctuations to the second power, vertical 
variations in the correlation timescale could still affect the nu- 
merical estimate of the diffusion coefficient. We thus calculated 
the value of Tcorr at two locations, one i n the disc m idplane and 
the oth er in the disc corona. Following iFromang & Papaloizoij 
(l2006l) . Tcorr Can be evaluated by monitoring the time variation 
of the function 



Szz(t) =< v,(z, t)v,(z, ^ + t) > , 



(23) 



where < . > denotes an ensemble average. Szz(t) is expected 
to decrease toward zero from initially positive values in a time 
Tcorr- To estimate the later, the model in which (Qt)o = 0.01 
was restarted at time = 573. We then calculated and averaged 
the function v^(z, to)Vz(z, to -\- r) at seven different radii R = 2, 
2.5, 3, 3.5, 4, 4.5 and 5 over a kernel of five cells in the radial 
direction. At all radii, the function was further averaged over 
two ranges in the meridional direction: \6\ < 1.5H/R to pro- 
duce the function S^^i, and \e\ e [1.5H/R, 3.5H/R] to produce 
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the function 8222 • ^zzl ^^^^ represents the velocity correlation 
function near the disc midplane while 8222 represents the corre- 
lation function in the upper layers including the corona. The two 
functions are plotted as a function of r in Fig. [S] respectively, 
with a solid and a dashed line (both curves are normalized by 
their value at r = 0). As expected, both display an initial de- 
crease toward zero around which they stabilize after a few tenths 
of an orbit^This qualitative trend is in agreement with the results 
of iFromang & Papaloizoul (l2006h . The dotted curves over plot- 
ted on the same figure represent the functions S^^ = exp(-^/To) 
for To = 0.05, 0.1, 0.15 and 0.2 orbits. Although it is difficult to 
measure the correlation timescale precisely given the large fluc- 
tuations we obtained, these curves can still be used as a way to 
estimate the correlation time Tcorr governing the functions S^-^i 
and 8222- They indicate that Tcorr is ~ 0.05 - 0.1 in the corona 
and ~ 0.1 - 0.2 in the disc midplane. B oth values are compara- 
ble to the value of 0. 15 orbit reported by Froman g & Papaloizoul 
(l2QQ6h . Although these two timescales are difl'erent, their ratio is 
at most two and certainly accounts for less variation in the dif- 
fusion coeflftcient D than the variations in the vertical velocity 
fluctuations described above. Therefore, when plugging numer- 
ical estimates for the correlation time into Eq. (|2T1) , we will use 
Tcorr = 0.15 orbits in the remaining of this paper. This value is 
such that Q^Tcorr ~ 1 and is i n agreement with prev iously re- 
ported results in the literature ( Johansen et al. '2006b') and with 
values used to model the efl'ect of turbul ence in theoreti c al stud - 
ies of planet formation (see for example lWeidenschillinglll984l) . 
thus giving additional support to such work. 

4.3. Dust spatial distribution 

For the three simulations we performed, we averaged the dust 
density in azimuth and over time between ^ = 550 and t = 600 
orbits. The spatial distributions we obtained using this procedure 
are shown in Fig. [6l for the models having (Qt)o = 0.01 (left 
panel), 0.001 (middle panel) and 0.0001 (right panel). Note that 
the radial extent of the snapshots is limited to R < 6. At larger 
radii, the outer bufl'er region we use (see Fromang & NelsonI 
[2006h starts to afl'ect the dust distribution. As expected, the 
smaller the dust particle sizes, the thicker the dust disc. This is 
simply because smaller particles are better coupled to the gas 
and can thus be lifted further away from the disc midplane by 
the turbulence before they decouple from the gas. The left panel 
also shows some sign of the dust disc flattening at large radii. 
This is simply because the strength of the turbulence, as mea- 
sured for example by the parameter a, decreases at large radii, 
as seen in Figs. [2] and [3l It should not b e confused with the ap- 
parent dust disc flattening identified by Dullemond & Dominik 
(120041) . which occurs because of self-shadowing in the presence 
of weak turbulence. Finally, Fig. [6l also highlights one of the lim- 
itations of our work: with the set-up of the simulations presented 
in this paper, we cannot easily extend our parameter survey to 
smaller particles. Indeed, when (Qt)o = 0.0001, the dust disc 
already covers almost the entire computational domain in the 
vertical direction. Reducing further the size of the dust particles 
would require an increase in the size of the computational do- 
main in the meridional region for the decoupling between gas 
and dust to occur within the computational grid. This would re- 
quire an increase in the computing time required for a simulation 
performed with the same resolution. This will soon become pos- 
sible as computational resources improve, but is currently very 
challenging. 

To compare these results with the models presented in 
Sects, im 13.2. II and 13.2.21 we computed the expected 2D dust 



distribution that each of them would predict. In doing so, we 
used the same disc and dust parameters as in the simulations. 

- For the first model (presented in Sect. 13.1b , we fitted a 
Gaussian profile to the vertical profile of the dust density 
at each radius. This was done by performing a least squares 
fit to the dust density profile over the entire vertical extent 
of the disc. The resulting spatial distributions are plotted in 
Fig. [7] using the same color table and spatial domain as in 
Fig.ia 

- For the second model (presented in Sect. 13.2.11) , we used 
Eq. ST9\i to calculate the vertical profile of the dust density 
at each radius. The numerical value of the dimensionless 
diff'usion coefficient D at each radius was calculated using 
Eq. ([2Qb in which we plugged the value of a calculated ac- 
cording to Eq. ([22b . We used Sc = 1.5 as this value turns out 
to provide the best fit to the simulations. Being of order unity, 
it is also in agreeme nt with the results of previous local nu- 
merical simulations ( Johansen & Kla hd2005l: l Johansen et al] 
2006b; Fro mang & Papaloizou 2006). The resulting spatial 
distributions of the dust density are plotted in Fig.[8l 

- For the third model (presented in Sect. 13.2.2b , we integrated 
Eq. ([TSb numerically, using at each disc altitude the estimate 
for the diff'usion coefficient D provided by Eq. ([2T1) . In this 
equation, we used the azimuthally and time averaged vertical 
profile of the vertical velocity fluctuations at each radius (as 
shown on Fig.lHfor the special case R = 2.93) and the cor- 
relation time Tcorr = 0.15 orbit as explained in Sect. 14.21 The 
resulting spatial distributions of the dust density are plotted 
in Fig. [3 

The spatial distributions of the dust density obtained with the 
three difl'erent models are in rough agreement with the simula- 
tions and with naive expectations: all predict that the dust com- 
ponent of the discs thicken as dust particles decrease in size, 
in agreement with the MHD simulations. For the model having 
(I1t)o = 0.01 (i.e. for the largest particles), the agreement be- 
tween the three models and the simulation is quite good. This 
can be understood easily: in this model, the dust scale-height is 
about 0.3 H. In other words, the solid particles concentrate close 
to the equatorial plane, where the turbulence properties are fairly 
homogeneous (the vertical velocity fluctuations only start to rise 
significantly above ~ 1.5//). Thus the diff'usion coefficient cal- 
culated according to Eq. ([2T1) is roughly constant and the sec- 
ond and third models (with constant and varying diff'usion co- 
efficient) give a similar result. Moreover, close to the equatorial 
plane (i.e. Z ^ //), the leading order expansion to Eq. ([T9b turns 
out to be Gaussian, which explains the similarities between the 
three models for large particles. 

On the other hand, models having smaller dust sizes, i.e. 
(Qt)o = 0.001 and 0.0001, show diff'erences between the diff'er- 
ent approaches. The Gaussian fit to the simulations always over- 
estimates the dust density in the disc corona (Z/H > 2), showing 
that in general the vertical profile of the dust density distribution 
is non Gaussian. On the other hand, the model having a constant 
diff'usion coefficient always underestimates the dust density in 
the disc corona. Only the model taking the vertical variation of 
the diff'usion coefficient into account gives a satisfactory fit to 
the simulations, especially in the disc upper layers. 

These diff'erences can be made more quantitative by compar- 
ing the vertical profile of the dust density between the simula- 
tions and the models. This is done in Fig.[TOl The left panel gath- 
ers the results corresponding to the case (^t)o = 0.01, the mid- 
dle panel shows results obtained when (flT)o = 0.001. Finally, 
the right panel shows results obtained when (flT)o = 0.0001. 
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Fig. 6. Dust density distribution in the (R,Z) plane in the numerical simulation for the model having (flT)o = 0.01 (left panel), 
(Qt)o = 0.001 (middle panel) and (Qt)o = 0.0001 (right panel). For all cases, the raw data have been first azimuthally averaged and 
then time averaged between ^ = 550 and t = 600. 
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Fig. 7. Dust density distribution in the (R, Z) plane computed by fitting a Gaussian vertical profile to the simulation data (azimuthally 
and time averaged) at each radii. As is the case for fig.[6l the left, middle and right panels respectively correspond to (Q^t)o = 0.01, 
0.001 and 0.0001. 




Fig. 8. Same as fig.|7]using the model that assumes D = constant. 



S.Fromang & R.Nelson: Dust settling in turbulent protoplanetary discs 



9 







Fig. 9. Same as fig.|7]using the model that assumes D = SvItc, 



In each panel, the solid curves plot the vertical profile of the 
dust density obtained in the MHD simulation by temporally, az- 
imuthally and radially averaging the results. The radial averag- 
ing is performed over the range R e [3,5]. The dot-dashed, 
dashed and dotted curves, respectively, correspond to the model 
using a Gaussian fit to the data, a constant diff'usion coeffi- 
cient and a vertically varying diff'usion coefficient. Again, the 
left panel demonstrates the good agreement between all ap- 
proaches in the case of large particles, as the curves used to rep- 
resent the results of the diff'erent models are almost undistigu- 
ishable on that plot. It also provides, however, a first hint that 
the Gaussian fit increasingly overestimates the dust density as 
one moves away from the midplane, even for these larger parti- 
cles. The middle and right panels confirm these results and show 
that the model having a vertically varying diff'usion coefficient 
provides in general the best fit to the data in the disc corona. 
For example, in the case (Qt)o = 0.001, the three models give 
a good fit to the data for Z < 2H. For Z > 2H, the Gaussian fit 
to the data overestimates the density obtained in the simulations 
(their ratio is about 10^ at Z = 3H) and the model having a con- 
stant diff'usion coefficient underestimates that density (the ratio 
at 3H is in excess of 10^). The agreement between the simula- 
tions and the model having a vertically varying diff'usion coef- 
ficient is better as the ratio between the predicted and observed 
dust densities is always bounded by 10, despite the fact that the 
dust density itself varies by more than 6 orders of magnitude. It 
is interesting to point out, though, that the value of the density 
predicted by this model close to the equatorial plane seems to 
significantly underestimate the results of the simulations. This is 
most likely because the correlation timescale we used in calcu- 
lating the diff'usion coefficient underestimates the midplane cor- 
relation timescale of the turbulence (see section 14.21) . Using a 
varying correlation timescale would certainly further improve 
the agreement with the numerical simulations, but such a level 
of refinement is probably meaningless given the other approx- 
imations involved in the simulations themselves. As shown on 
Fig.[TOl the situation is similar in the case (^t)o = 0.0001: the 
best fit to the numerical simulations is provided by the final and 
more elaborate model. Indeed, at Z > 2.5//, the dust density pre- 
dicted by the model with a constant diff'usion coefficient starts to 
underestimate significantly the results of the numerical simula- 
tions. Above Z = 3//, the diff'erence becomes enormous, as in 
the case (^t)o = 0.001. The agreement, however, seems better 



with the model that uses a Gaussian fit to the data than in the 
case (Qt)o = 0.001. The Gaussian fit indeed gives a good fit 
up to Z = 3H. This apparent agreement would most probably 
break down for Z > 4.3, as the ratio with the simulated densities 
is already greater than two orders of magnitude at Z = 4H and 
increases with increasing Z. 

5. Discussion and conclusions 

In this paper, we have studied dust settling in turbulent proto- 
planetary discs using global MHD numerical simulations per- 
formed with the code GLOBAL using spherical coordinates. In 
this section, we summarize our main findings and discuss the 
limitations and future prospects of our work. 



5.1. Dust density vertical profile 

The first point that emerges is that Gaussian profiles fail badly 
to reproduce the data extracted from the simulations. This point 
is further illustrated by Fig. [TT] where we compare the dust pro- 
file in the case (^t)o = 0.001 (solid line) with a set of Gaussian 
profiles (dotted lines) computed according to the following equa- 
tion. 



f 



Pd,th = exp 



2/^L. 

d,th 



(24) 



in which we used the values Hd,th/H = 0.4, 0.5, 0.6, 0.7, 0.8 and 
0.9. None of the dotted lines is an acceptable fit to the data. Close 
to the disc midplane, by eye inspection seems to suggest that the 
dust disc scale-height is ~ 0.9// while one would estimate it to 
be of order 0.6// at Z ~ 3H. Both values are diff'erent from the 
dust disc scale-height returned by the least squares fit to the data, 
H^/jH = 0.84. The physical reason for this behaviour is clear: 
close to the midplane, where the gas density is high, the gas- 
dust coupling is strong and the dust traces closely the Gaussian 
profile of the gas, whereas in the disc upper layers the coupling 
is weak and the dust-gas ratio decreases as one moves further 
away from the midplane. These results emphasize the point that 
estimates of the dust disc scale-height obtained using a Gaussian 
fit may lead to incorrect conclusions. 

Nevertheless, in order to compare with results published pre- 
viously in the literature, we measured the dust disc scale-height 
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Fig. 10. Vertical profiles of the dust density, radially averaged 
between R = 3 and R = 5. As shown by the legend displayed 
on the upper panel, the diff'erent curves (normalized by their 
midplane values) were obtained using the numerical simulations 
(solid line), by fitting a Gaussian profile to the simulation data 
(dotted-dashed line), using a constant turbulent diff'usion coeffi- 
cient (dashed line) and a vertically varying diff'usion coeflftcient 
(dotted line). The upper, middle and bottom panels respectively 
correspond to (Qt)o = 0.01, 0.001 and 0.0001. In the first case, 
all models successfully reproduced the simulation, while for the 
smaller dust sizes, only the later and more elaborate model pro- 
duces satisfactory results. 




Fig. 11. This figure compares the dust density vertical profile 
(shown with a solid line), in the case (Q^t)o = 0.001, with a set of 
Gaussian profiles with diff'erent scale-heights Hdjh (shown with 
dotted lines). The diff'erent values of Hd,th we used are 0.4//, 
0.5H, 0.6H, O.IH, 0.8// and 0.9H. None of the dotted lines 
provides an acceptable fit to the numerical dust density. This 
demonstrates that the dust vertical profile is not Gaussian. 
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Fig. 12. Dust disc scale-heights obtained through a Gaussian fit 
to the data are represented by the diamonds as a function of 
(Qt)o. The solid line displays the best power law fit to the points 
and indicates that Hd/H oc a~^-^. 



in our simulations using such a fit. We found Hd,th/H = 1.0, 
0.84 and 0.40, respectively, for (Dt)o = 0.0001, 0.001 and 0.01. 
The variation of Hd,th/H as a function of (^t)o is shown by the 
diamonds in Fig. [121 The dotted line displays the function 



Hd/H = 0.7 



(^T)o 

0.001 



-0.2 



(25) 



which is the best fit to the data. Therefore, if we were to anal- 
yse our simulations assuming that the dust density is Gaussian, 
we would obtain oc a~^-^ since (Qt) q qc a. Inter e stingly , this 
is not too diff'erent from the results of iPinte et al.l (l2008h who 
report an exponent equal to -0.05, although it is clear that our 
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results show a stronger relationship between the dust disc scale- 
height and particle size. Both values are, however, largely differ- 
ent from the value -0.5 which is often reported in the literature 
(iDubruUe et al.[l995l:ICarballido et al.l2QQ6D . This is because the 
latter is obtained when solving Eq. (fTTl) in the strong settling 
limit Hd ^ H that is mostly relevant for large particles. In this 
case, the vertical variation of Qr^ can be neglected, the vertical 
profile is Gaussian and the exponent -0.5 is recovered. For the 
small particles we study here, the dust disc is thick and these ver- 
tical variations have to be taken into account, which leads to the 
more complicated expression given by Eq. ST9\j and departure 
from a Gaussian profile. 

5.2. The Schmidt number 

The second result that emerges from our work is that a model 
having a constant diff'usion coefficient increasingly underesti- 
mates the dust density as the particle size a is decreased. In 
other words, the vertically averaged dust diff'usion coefficient 
decreases with a. This is not unexpected, since the Schmidt 
number introdu ced in Sect. 13.2. II is known to be an increasing 
function of Qr (ICiizzi et al.lll993l: [S chrapler & Henning ll2004t 
lYoudin & Lithwickll2007l) . It is unclear, however, whether such 
studies, which assume homogeneous Kolmogorov-like turbu- 
lence, are applicable to the highly magnetised flow of the corona. 
This is why we did not attempt to make a direct comparison be- 
tween our results and these theories. For the sake of complete- 
ness, it is nevertheless instructive to report here the vertically av- 
eraged Schmidt number we measured in each case. As described 
above, Sc = 1.5 already provides a good fit of the dust density in 
the case (Qt)o = 10"^. For the cases (I1t)o = 10"^ and 10""^, we 
found that the vertically averaged Schmidt numbers that best fit 
the data are respectively Sc = 0.4 and Sc = 0.030. Although they 
indicate a scaling with (flT)o close to linear, t hese results are not 
necess ary in disagreement with the result of Youdin & Lithwickl 
iZOOl), who found a quadratic scaling, because of the vertical 
average we made when doing such measurements. Note also that 
these fairly low values of the Schmidt number indicate that tur- 
bulent diff'usion of the smallest dust particles can be much more 
efficient than angular momentum transport. This diff'erence orig- 
inates from their diff'erent physical origin: angular momentum 
is transported radially by the Maxwell and Reynolds stresses 
while dust diff'uses away from the disc midplane because of 
the vertical velocity fluctuations of the gas. The former quan- 
tities decrease as one moves away from the dis c midplane, (see 
iMifler & Stondl2000l: 'Froman g & Nels'onll2006h while the latter 
quatity increases away from the disc midplane. This means that 
angular momentum is inefficiently transported in the disc corona 
while small dust particles are efficiently diff'used at the same lo- 
cation. As a consequence, the Schmidt number (i.e. the ratio of 
both diff'usion coefficients) is much smaller than unity for the 
smallest of our particles. 




Fig. 13. Comparison between the dust density vertical profile 
(shown in both panels with the solid line) in the case (0^t)o = 
0.001 with a simple toy model in which the turbulent velocity 
distribution is a step function (see text for details). In the upper 
panel, the velocity fluctuations are taken to be dv^^upl^s -0.15 
for \z\ > 2H, while Svz,mid/<^s = 0.025 (dotted line), 0.05 (dashed 
line) and 0.075 (dotted-dashed line) for |z| < 2H. In the lower 
panel, the velocity fluctuations are Svz,mid/(^s = 0-05 for \z\ < 2H 
and Svz,up/Cs ~ 0.075 (dotted line), 0.15 (dashed line) and 0.30 
(dotted-dashed line) for \z\ > 2H. The results show that the dust 
density vertical profile is fairly insensitive to the midplane ve- 
locity fluctuations but more strongly depends on their amplitude 
in the disc upper layers. 



5.3. A toy model 

The main result of this paper is the construction of a simple 
model that gives a reasonable fit to the simulations. In this 
model, the dust particles are diff'used away by turbulence with 



^ The fact that we obt ain S c hmidt numbers lowe r than one, 
in contrast to Cuzzi etjO] (Il993h . ISchrapler & Hennind (2004) and 
00% i 

to that used in these theoretical studies 



lYoudin & Lithwick (2007 



lYoudin & Lithwickl (120071) 



is due to our definition being different 
as discussed in detail by 



a diff'usion coefficient that scales with the square of the turbulent 
velocity fluctuations. Accordingly, it should be possible in prin- 
ciple to extract the vertical profile of the velocity fluctuations 
from the dust density vertical profile. This is, however, an inver- 
sion problem. As such, it is susceptible to being degenerate and 
we shall see in this section that this is indeed the case. 

The relevant question to ask in this context is the following: 
provided we are able to measure the dust density vertical pro- 
file, what is it mostly sensitive to? Can we hope to constrain the 
velocity fluctuations in the disc midplane or is it mostly a conse- 
quence of the amplitude of the fluctuations in the upper layers. 
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To answer that question, we designed the following toy model: 
guided by the vertical profile of the velocity fluctuation ampli- 
tudes shown in Fig.jH we considered the following analytic ver- 
tical profile for the velocity fluctuations: 

^ I ^Vz,mid + V^Vz,up - ^Vz,mid] {^f if kl < 2H ^26) 

1 Svz,up Otherwise 

where Sv^^mid and Svz,up stand for the turbulent velocity fluctu- 
ations in the disc midplane and in the disc corona. As shown 
on Fig. m typical numerical values are Svz,mid/Cs ~ 0.05 and 
^^z,up/Cs ~ 0.15. To investigate the sensitivity of the results to 
the midplane velocity fluctuations, we calculated the dust den- 
sity vertical profile in the case (^t)o = 0.001 by numerically 
integrating Eq. ([TSl) , using Eq. (|2T1) and (l26b with Svz,mid/Cs = 
0.025, 0.05 and 0.075 and Svz,up/Cs = 0.15 (i.e. we used in the 
disc corona the value suggested by the simulation data) . The re- 
sults are summarized on the upper panel of Fig. \T3\ where the 
dust density profile in the case (^t)o = 0.001 is shown with 
the solid line, while the numerically integrated profiles are rep- 
resented by dotted, dashed and dot-dashed lines for Svz,mid/(^s = 
0.025, 0.05 and 0.1, respectively. The last three curves are very 
similar in this plot and all give a fairly good fit to the simulations, 
especially in the disc upper layers. This shows that the dust den- 
sity vertical profile is fairly insensitive to the midplane velocity 
fluctuations. To estimate the sensitivity of the dust density pro- 
file to the velocity fluctuations in the upper layers, we repeated 
the same analysis using Svz,mid/(^s = 0.05 (i.e. we used in the 
disc midplane the value suggested by the simulation results) and 
^^z,up/(^s = 0.075, 0.15 and 0.30. The results are shown on the 
bottom panel of Fig. [T3l with the same conventions as for the up- 
per panel. In this case, the dust density vertical profile is seen to 
be much more sensitive to the upper layer velocity fluctuations 
and only the dashed curve, for which dv^^uplcs = 0.15 (i.e. in 
rough agreement with the numerical data), is in good agreement 
with the data. 

The implications of these results are twofold. First, the sim- 
ple toy model described by Eq. (l26l) would be suitable to use 
when trying to fit the observations as it reproduces the numer- 
ical data fairly well if we choose the values dv^^midlcs = 0.05 
and Svz,up/Cs = 0.15, compatible with the simulations. But these 
results also show that such a fit would only provide sensitive 
information about the turbulent velocity fluctuations in the disc 
upper layers. The physical reason for this last point is that such 
a fit is mostly sensitive to the properties of the region where the 
gas and dust decouple. For the small particles studied in this pa- 
per (and observed using the Spitzer telescope), this region turns 
out to lie in the disc upper layers. When observing larger par- 
ticles at longer wavelengths (for example with ALMA), it will 
become possible to constrain the turbulent velocity fluctuations 
of the disc closer to the midplane as such particles will settle and 
decouple deeper in the disc. 

5.4. Limitations and future prospects 

Of course, there are strong limitations to our work due to the 
complex and CPU-intensive nature of global simulations of pro- 
toplanetary discs. On the purely numerical side, the limited res- 
olution we used is of course an issue, as was pointed recently in 
a nu mber of studies (iFromang & Papaloizoull2007l: ISimon et aP 
'2008). Proper simulations should include explicitly microscopic 
difl'usion coefficients (viscosity and resistivity), as the latter have 
been shown to be imp ortant in determining the saturation level o f 
the turbulence (,Fromang et al.ll200"7l: iLesur & Longarettill2007li) . 



However, the resolutions required to include these processes in 
global simulations are currently out of reach and one must in- 
stead rely on the subgrid model provided by numerical dissi- 
pation to carry out global numerical simulations. On a more 
physical side, there are also limitations due to the simple disc 
model we used. The locally isothermal equation of state we used 
is not appropriate for the disc inner parts that we are simulat- 
ing as the gas there is o ptically thick. Thi s co uld have numer- 
ous eff'ects. For example jPuUemondl (l2002h and iD' Alessio et al.] 
(1998) report a temperature increase in the inner disc upper lay- 
ers. Such an increase would strengthen the coupling between gas 
and grains in the disc corona (through a decrease in the local 
value of the parameter Qr^). This would cause small particles to 
settle less than reported in this paper and could change the re- 
lationship between and a. Obviously, the significance of the 
comparison we te ntatively made bet ween our simulations and 
the observations of iPinte et al.l (|2008) should be taken with care. 
Another topic of concern in our simulations is the assumption of 
ideal MHD. It is indeed well known that protoplanetary discs are 
so poorly ionized because of their large densities and low tem- 
peratures that parts of the flow, refered to as dead zones, remain 
laminar ( Gammia 1996). We completely ignored the efl'ects of 
dead zones in the present paper. Clearly, future work should im- 
prove the thermodynamic treatment of the gas, possibly includ- 
ing radiative transfer and dead zones. 

Nevertheless, we have shown in this paper that dust obser- 
vations can be used in principle to constrain the properties of 
MHD turbulence in discs. We have found that even the simplest 
simulations provide disagreements with previously used fits and 
diff'usion models because of the nature of disc turbulence. This 
illustrates even further the need to compare directly observations 
and numerical simulations. It will be important in the future to 
generate a grid of more realistic discs models (varying the disc 
parameters, including dead zones, flaring discs, non isothermal 
discs) and produce synthetic observations that could be com- 
pared in the next few years with multiwavelengths observations. 
This comparison would provide diagnostics of disc turbulence, 
the existence (or not) of dead zones and thus constrain planet for- 
mation models. The recent work of lPinte et al.l (l2008h shows that 
such multiwavelengths observations are starting to become fea- 
sible. When combined with future instruments like Herschel and 
ALMA, large samples will become available and will provide 
a wealth of constraints on disc structure and properties when 
combined with appropriate global numerical simulations of pro- 
toplanetary discs. 
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